# In this script we perform:

#  1. full asset-liability matching using 15 bonds, with maturities spanning over
#     each of the 15 time periods - the ones used in script
#     "S3.Estimating the market yield curve by replication.R".

#  2. cash flow immunization using 3 bonds, defined below.

library(tidyverse)
library(readxl)

############### Example of Asset-Liability Matching #############

# Load the market yield curve and the matrix of payments of bonds, taken from
# the same script.

load("Examples/market.yield.curve.Rdata")
market.yield.curve
payments.tbl <- read_xlsx(
  "Examples/interestRateData2026.xlsx", sheet = "paymentMatrixC"
)

# Consider matrix C
C <- as.matrix(payments.tbl)
C

# Assume a liability cash flow of 1M per time period.
data <- market.yield.curve %>% mutate(liab.cf = 1000000)
data

# Let us find the weights of the (full asset-liability) matching portfolio:
# (w_1,...,w_n) = (L_1,...,L_n) * inverse(C)
w.match <- data$liab.cf %*% solve(C) %>% as.numeric
w.match

# Verification
w.match %*% C

#################### Example of Cash Flow Immunization ################

# Most times, cash flow matching is not practical.
# We will perform cash flow immunization using 3 bonds.
maturities <- c(2, 8, 11)
coupons <- c(0.04, 0.05, 0.05)
face.values <- c(100, 100, 100)

# Let us add the bond payments to our data
data <- data %>%
  # Include the coupon payments
  mutate(
    coupons.1 = if_else(t <= maturities[1], coupons[1] * face.values[1], 0),
    coupons.2 = if_else(t <= maturities[2], coupons[2] * face.values[2], 0),
    coupons.3 = if_else(t <= maturities[3], coupons[3] * face.values[3], 0)
  ) %>% 
  # Include the principal payments
  mutate(
    principal.1 = if_else(t == maturities[1], face.values[1], 0),
    principal.2 = if_else(t == maturities[2], face.values[2], 0),
    principal.3 = if_else(t == maturities[3], face.values[3], 0)
  ) %>%
  # Compute total payments
  mutate(
    payments.1 = coupons.1 + principal.1,
    payments.2 = coupons.2 + principal.2,
    payments.3 = coupons.3 + principal.3
  ) %>% 
  select(-starts_with("coupons"), -starts_with("principal"))

# select(-c(coupons.1, coupons.2, coupons.3))

# Compute the duration and convexity of the liabilities cash flow (in the same 
# way as in script "4.Estimating Present Value Sentitivity to Yield Changes.R")
sens.liab <- data %>% 
  mutate(PV = P_t * liab.cf) %>% 
  mutate(
    d_aux = t * P_t * liab.cf,
    c_aux = t^2 * P_t * liab.cf
  ) %>% 
  summarize(
    PV = sum(PV),
    d_aux = sum(d_aux),
    c_aux = sum(c_aux)
  ) %>% 
  mutate(duration = d_aux/PV, convexity = c_aux/PV) %>% 
  select(PV, duration, convexity)
sens.liab

# Now we compute the duration and convexity of each of the 3 bonds
# We will encapsulate our code in a function in order to avoid repetition when 
# computing duration and convexity of our three bonds. Please note that when
# we invoke a column from a data frame inside a function using tidyverse, we
# have to "embrace" it, meaning that we surround it with {{...}}.

f.sensitivity <- function(table, variable){
  
  table %>% 
    mutate(PV = P_t * {{variable}}) %>% 
    mutate(
      d_aux = t * P_t * {{variable}},
      c_aux = t^2 * P_t * {{variable}}
    ) %>% 
    summarize(
      PV = sum(PV),
      d_aux = sum(d_aux),
      c_aux = sum(c_aux)
    ) %>% 
    mutate(duration = d_aux/PV, convexity = c_aux/PV) %>% 
    select(PV, duration, convexity)
  
}

sens.bond.1 <- f.sensitivity(table = data, variable = payments.1)
sens.bond.2 <- f.sensitivity(data, payments.2)
sens.bond.3 <- f.sensitivity(data, payments.3)
sens.bond.1
sens.bond.2
sens.bond.3

# In order to find the correct amounts of each bond to include in our portfolio W,
# we will construct the system of equations prescribed for second order immunization:
# PV_L = PV_A                             (equate present values)
# PV_L * Duration_L = PV_A * Duration_A   (equate dollar durations)
# PV_L * Convexity_L = PV_A * Convexity_A (equate dollar convexities)

# Function for computing Dollar Duration and Dollar Convexity
f.dollar <- function(table){
  
  table %>%
    mutate(DD = PV * duration, DC = PV * convexity) %>% 
    select(PV, DD, DC)
  
}

PV.sens.liab <- f.dollar(sens.liab)
PV.sens.liab
PV.sens.bond.1 <- f.dollar(sens.bond.1)
PV.sens.bond.1
PV.sens.bond.2 <- f.dollar(sens.bond.2)
PV.sens.bond.2
PV.sens.bond.3 <- f.dollar(sens.bond.3)
PV.sens.bond.3

# Our system of equations is of the form M * w = v, where:
#  - M is a matrix corresponding to the PV and sensitivity of the assets (known)
#  - w is the vector of portfolio weights (unknown)
#  - v is a vector corresponding to the PV and sensitivity of the liabilities (known)

# t() is the 'transpose' function in base R.
t(PV.sens.bond.1)
t(PV.sens.bond.2)
t(PV.sens.bond.3)
M <- cbind(t(PV.sens.bond.1), t(PV.sens.bond.2), t(PV.sens.bond.3))
M
v <- t(PV.sens.liab)
v

# Solving the system of equations (w = inverse(M) * v):
w <- solve(M) %*% v
w

# These weights are our solution. We have determined an immunizing portfolio to
# our liabilities. Let us check the asset cash flow, as well as its sensitivity
# measures.
data
data <- data %>% 
  mutate(
    # Cash flows for each bond
    cf.1 = payments.1 * w[1],
    cf.2 = payments.2 * w[2],
    cf.3 = payments.3 * w[3],
    # Assets cash flow
    assets.cf = cf.1 + cf.2 + cf.3
  )
# Same D and C as liabilities? Please check!
data %>% f.sensitivity(assets.cf)

# Same DD and DC as liabilities? Please check!
data %>% f.sensitivity(assets.cf) %>% f.dollar() 

# Calculate the present value of the surplus
data <- data %>%
  mutate(
    surplus =  assets.cf - liab.cf,
    PV.liab = P_t * liab.cf,
    PV.assets = P_t * assets.cf,
    PV.surplus = P_t * surplus
  )

data %>% summarize(sum(PV.surplus))

# Save the immunization data for later
immuniz.data <- data %>% 
  select(t, liab.cf, assets.cf, surplus, starts_with("PV"))
save(immuniz.data, file = "Examples/immuniz.data.Rdata")

totals <- data %>%
  summarize(
    across(c(liab.cf, assets.cf, surplus, PV.liab, PV.assets, PV.surplus), sum)
    )
# The PV of our surplus is approximately zero
totals

# Finally, we will see how our portfolio holds if we add a few perturbations
# to our yield curve. We include a perturbation 'a' for the first three years
# and a perturbation 'b' for the remaining years.
a <- 0.005
b <- -0.005
yield.perturbation <- c(rep(a,3),rep(b,9))
yield.perturbation

# Homework: compute the present value of the liabilities, assets and surplus in
# this new case. It is easy!

perturbed.data <- data %>% 
  select(t, y_t = y_t.cont, liab.cf, assets.cf) %>% 
  mutate(
    y_t = y_t + yield.perturbation,
    P_t = exp(-t * y_t)
    ) %>% 
  mutate(
    liab.PV = P_t * liab.cf,
    assets.PV = P_t * assets.cf
  )

# PV of assets, PV of liabilities and PV of surplus
perturbed.data %>%
  summarize(liab.PV = sum(liab.PV), assets.PV = sum(assets.PV)) %>% 
  mutate(surplus.PV = assets.PV - liab.PV)






# Conclusion: this method provides protection against "parallel" shifts in the
# yield curve, but not so much against "mixed" perturbations.